Mooney–Rivlin Parameter Determination Model as a Function of Temperature in Vulcanized Rubber Based on Molecular Dynamics Simulations

The automotive industry is entering a digital revolution, driven by the need to develop new products in less time that are high-quality and environmentally friendly. A proper manufacturing process influences the performance of the door grommet during its lifetime. In this work, uniaxial tensile tests based on molecular dynamics simulations have been performed on an ethylene–propylene–diene monomer (EPDM) material to investigate the effect of the crosslink density and its variation with temperature. The Mooney–Rivlin (MR) model is used to fit the results of molecular dynamics (MD) simulations in this paper and an exponential-type model is proposed to calculate the parameters C1(T) and C2T. The experimental results, confirmed by hardness tests of the cured part according to ASTM 1415-88, show that the free volume fraction and the crosslink density have a significant effect on the stiffness of the EPDM material in a deformed state. The results of molecular dynamics superposition on the MR model agree reasonably well with the macroscopically observed mechanical behavior and tensile stress of the EPDM at the molecular level. This work allows the accurate characterization of the stress–strain behavior of rubber-like materials subjected to deformation and can provide valuable information for their widespread application in the injection molding industry.


Introduction
Elastomers are widely used in industry because of their exceptional properties.These include elasticity and durability.Synthetic rubber components are essential in automotive parts for sealing systems.Typically, elastomers are used in applications where they deform at low strain rates.Ethylene-propylene-diene monomer (EPDM) rubber is chemically stable; this type of elastomeric polymer is characterized by high elasticity and resistance to chemical media, while the fatigue strength and thus the service life of EPDM components are affected by increasing temperature and load cycles [1].
The automotive industry is entering a digital revolution, driven by the need to rapidly develop new products that integrate more efficiently and sustainable automated and autonomous technologies.In recent years, the industry involved in the production of EPDM components has implemented sustainability and carbon footprint reduction policies; to achieve these objectives, sustainable practices have been implemented for the design and manufacture of these components, using virtual prototyping tests to develop quality and environmentally friendly products [2].
The effectiveness of the solutions resulting from these numerical methods depends on the accuracy with which the mechanical properties of the rubber material can be explained, so numerous theoretical models have been developed in an attempt to describe the mechanical behavior of rubber.The macromolecular network structure of elastomers allows these materials to undergo large deformations when subjected to tensile stress, resulting in a nonlinear relationship between stress and strain.The effect of crosslinking and the temperature during the curing process is considered one of the major challenges in studying the complex mechanical behavior of amorphous crosslinked materials [3][4][5].
To better understand the influence of temperature on crosslinking and the stress-strain relationship, several physical explanations have been summarized, ranging from chain breakage in the amorphous structure to the sliding of molecules and the formation of more complex composite structures.For the complex crosslinking process during curing, it is evident that no existing model can accurately and effectively capture the dependence of the crosslink density and temperature on the characteristic tensile behavior of an elastomer in a thermodynamically consistent and numerically robust manner due to the strong nonlinearity [6].
Current models for the stress-strain (S-s) relationship consider only the idealized phenomenon, without taking into account intrinsic phenomena such as the crosslink distribution or crosslink density.Researchers in this field have proposed numerous hyperplastic constitutive models to describe the mechanical behavior of elastomers under stress; these models can be divided into two categories: The first category is based on statistical mechanics and assumes that the material consists of randomly distributed molecular chains.The second category is the phenomenological model, which applies the continuum mechanics method and assumes that the elastic potential is formed by invariants constructed by the strain or S-s ratio [7].
The properties of viscoelastic materials can be characterized as a function of the strain energy density.Models have been developed to theoretically reproduce S-s curves from experimental strain data.The mechanical behavior exhibits a high degree of nonlinearity; therefore, it is complex to establish a generalized strain energy density function that competently describes the S-s relationship experimentally.
The crosslink density is an important parameter affecting the mechanical properties of elastomeric materials, and these properties change significantly with the increasing crosslink concentration.Thermal gradients during the rubber vulcanization process strongly influence the density and distribution of crosslinks and, therefore, have a decisive influence on the mechanical properties and the final performance of the material [8].The viscosity and elasticity of the crosslinked polymers are strongly dependent on the morphology of the structure and the length of the chain segments, as this determines the dynamics.Theoretically, adjacent chains in the amorphous structure can only move locally with the constraints imposed by the network bonds, leading to a reduction in the configurational freedom of mobility as well as an increase in local friction and intermolecular coordination.
The elasticity of rubber is mainly determined by the formation of chemical and physical crosslinks during the curing process, where the chain length and the crosslink density define the dynamics and the stored elastic entropy [9].
In phenomenological models, it is assumed that the elastic potential is formed by invariants constructed by the strain ratio or strain, and the different parameters of the elastic potential are determined based on experimental data using, for example, the Ogden model, the Yeoh model, or the Mooney-Rivlin model.In the Mooney-Rivlin model, the crosslink density is an important parameter that affects the S-s ratio; the parameters in this model are associated with intrinsic material properties such as network structures and intermolecular forces [10][11][12].
Several investigations have been reported using a molecular dynamics simulation on polymer properties, which explore the kinetic properties of the EPDM considering the free volume and crosslink density.For example, Lengger et al. [11] use the Mooney-Rivlin model with photoelasticity for experimental measurements to determine the evolution of cure-dependent material parameters.Vishvanatperumal et al. [13] determine the crosslinking density by studying the experimental stress-strain behavior with finite element analysis by fitting the experimental data to the Mooney-Rivlin model.Morovati et al. [14] studied the effect of fatigue stress on elastomers, providing information on the complexity of the constitutive behavior of the crosslinked polymers.Liu et al. [12] proposed a hyperplastic constitutive model that shows good agreement with the experimental data obtained in elastomeric materials subjected to large deformations.Using molecular dynamics simulations, Wang Y. et al. [5] found that the crosslinking produced in the side chains influences the stiffness and enhances the inhibitory effect on the diffusion properties; the results of the study allow for an advance in the understanding of the microscopic aspects underlying the performance of EPDMs.Wang Ao et al. [15] used time-temperature superposition and varying strain rates in molecular dynamics simulations and the data obtained from the simulations were used to predict the constants of Rouse's exponential phenomenological model to determine the rheological properties of EPDMs.
It is important to note that the methodology of constructing crosslinked models using MD simulations has undergone rapid development in recent years.Initially, Yarovsky and Evans [16] proposed a static approach in 2002, in which the crosslinked bonds were created at a time forming a low-molecular-weight crosslinked network based on a predefined fixed radius.In this context and years later, Wu and Xu [17] created a small-scale crosslinked model using a dynamic method of iteratively creating new bonds with a fixed shear reaction radius by performing energy minimization and molecular dynamics (MM/MD) after each reaction.On the other hand, Varshey et al. [18] simulated the crosslinking reaction by proposing a stepwise method by constantly increasing the shear reaction radius, followed by MM/MD simulations, forming a relatively large crosslinked model.In general, to simulate the crosslinking process of EPDM, crosslinking should be added periodically by implementing some stepwise or multi-step algorithm as this is key to simulate the crosslinking process of EPDMs; the most prominent works in this regard to mention a few have been Papanikolaou et al., Wang et al. and van Duin et al. [2,19,20].
In this research, the study of the S-s relationship using the Mooney-Rivlin model considering the crosslink density in the amorphous structure of EPDMs through molecular dynamics modeling is approached.The methodology proposed in this work considers the molecular dynamics of the polymer chains in the amorphous structure on the macroscopic behavior for the mechanical characterization of EPDM subjected to large deformations and to provide a phenomenological model for their design in applications for the automotive industry.

Molecular Dynamics Model
Molecular dynamics (MD) simulations assume a non-overlapping configuration of all amorphous EPDM polymer chains, with five chains consisting of 1000 monomers randomly distributed in a simulation box (Figure 1a).Three-dimensional periodic boundary conditions are applied to eliminate surface effects.
The force field implemented in this work was calculated by condensed-phase optimized molecular potentials for atomistic simulation studies COMPASS.A simplified model of the Coulomb and van der Waals force interactions is used for the COMPASS force field [21,22].The van der Waals interactions hold atoms of opposite charge at a certain distance from each other.The summation methods for electrostatic (Coulombic) and van der Waals interactions were calculated with the Atom-based and Ewald algorithms, respectively [19,21].
The velocity-Verlet algorithm is employed to integrate the equations of motion, with a time step integration set to 1 fs, and the initial velocities were chosen randomly from the Boltzmann distribution.An Andersen-type thermostat and barostat were used to control the temperature and pressure.
There is no general methodology for equilibrating the MD systems, but there are guidelines to follow (see Figure 1b).In this study, using the methodology proposed by Gómez et al. [4], the system was equilibrated in two stages with four frames (replicas), optimizing the topology for each frame (pre-equilibration).
The generated frames are equilibrated through various conditions to bring them to a constant energetic state, defining conditions related to variables such as energy, pressure, volume, and temperature.
In the first equilibrium stage, simulated annealing with a constant energy and volume NVE ensemble is carried out in a temperature range of 350 to 500 K with 100 cycles and 10 ramps per cycle, a time step 1 fs, with 200,000 steps.
In the second step, the system was relaxed at a constant volume and controlled temperature NVT ensemble using the Andersen thermostat.On the obtained trajectories, to equilibrate the pressure variation, NPT ensemble simulations (100 MPa) were performed to equilibrate the simulation box size and to obtain the correct system density for production, both simulations with 20 ns with a time step of 1 fs, and the process was repeated in the temperature range from 428 K to 478 K with 10 K steps.
The method described below is based on earlier work by Papanikolaou et al. and adopted by Wang Y et al. [2,5].The simulation was based on the calculation of all crosslinks and the distances between the reactive atoms within a preset cut-off range between 0.4 nm and 1 nm and a crosslinking limit of 90%.The crosslinking kinetics were performed randomly and periodically and covalent terms were created for crosslinking.The new bonds were activated by alternating equilibrium distances and force coefficients to avoid missing atoms in the simulation box.The lower bound cut-off distance was increased by 0.025 nm at each iteration, and the distance between the reactive atoms was rechecked.The simulation box was equilibrated using three-dimensional periodic boundary conditions to eliminate surface effects on dynamic crosslinking.The procedure continues until the reactive The interaction potential function in a system between polymer chains and its parameters associated with the valence terms represent the internal bond coordinates (b), angle (θ), torsion angle (φ), and out-of-plane angle (χ).The vibrational frequencies and the structural variations associated with the changes in their conformation are related to each other in terms of cross-coupling.These terms include combinations of two or three internal coordinates considering the bond angle potential (E ba ), the in-plane dihedral rotation potential (E dh ), the out-of-plane vibrational potential (E inv ), and the cross term (E cross ).Non-bonding interactions are described by Lennard-Jones 9-6 potential for the van der Walls term (E vdW ) and a Coulomb function for an electrostatic interaction (E coul ).The functional terms of COMPASS are given by Equation ( 1): where E total is the total energy of the system and E vdW is the potential energy due to intermolecular forces calculated by Equation ( 2): where the subscripts ij in the summations represent the potential energy interactions, ε ij represents the potential energy parameter, σ ij represents the finite distance when the potential energy between the molecules becomes zero, r ij denotes the cut-off distance between atoms i and j and E coul is the Coulomb electrostatic potential that can be expressed by Equation (3): where q i and q j are the fixed partial charges between atom i and j within the same molecule.E ba is the bond angle potential shown in Equation ( 4): where θ represents the energy associated with the angle, the value θ 0 represents the equilibrium, and k, k 1 , k 2 , and k 3 are force constants.E dh is the dihedral rotation potential within the molecule defined by Equation ( 5): where ϕ represents the torsional angle.The out-of-plane vibrational potential E inv is represented by Equation ( 6) Wilson out-of-plane χ internal coordinates are shown in Equation ( 7) E cross , which includes the cross-coupling terms to calculate the vibrational frequencies and structural variations associated with conformational changes.
The velocity-Verlet algorithm is employed to integrate the equations of motion, with a time step integration set to 1 fs, and the initial velocities were chosen randomly from the Boltzmann distribution.An Andersen-type thermostat and barostat were used to control the temperature and pressure.
There is no general methodology for equilibrating the MD systems, but there are guidelines to follow (see Figure 1b).
In this study, using the methodology proposed by Gómez et al. [4], the system was equilibrated in two stages with four frames (replicas), optimizing the topology for each frame (pre-equilibration).
The generated frames are equilibrated through various conditions to bring them to a constant energetic state, defining conditions related to variables such as energy, pressure, volume, and temperature.
In the first equilibrium stage, simulated annealing with a constant energy and volume NVE ensemble is carried out in a temperature range of 350 to 500 K with 100 cycles and 10 ramps per cycle, a time step 1 fs, with 200,000 steps.
In the second step, the system was relaxed at a constant volume and controlled temperature NVT ensemble using the Andersen thermostat.On the obtained trajectories, to equilibrate the pressure variation, NPT ensemble simulations (100 MPa) were performed to equilibrate the simulation box size and to obtain the correct system density for production, both simulations with 20 ns with a time step of 1 fs, and the process was repeated in the temperature range from 428 K to 478 K with 10 K steps.
The method described below is based on earlier work by Papanikolaou et al. and adopted by Wang Y et al. [2,5].The simulation was based on the calculation of all crosslinks and the distances between the reactive atoms within a preset cut-off range between 0.4 nm and 1 nm and a crosslinking limit of 90%.The crosslinking kinetics were performed randomly and periodically and covalent terms were created for crosslinking.The new bonds were activated by alternating equilibrium distances and force coefficients to avoid missing atoms in the simulation box.The lower bound cut-off distance was increased by 0.025 nm at each iteration, and the distance between the reactive atoms was rechecked.The simulation box was equilibrated using three-dimensional periodic boundary conditions to eliminate surface effects on dynamic crosslinking.The procedure continues until the reactive pairs are exhausted within the pre-assigned cut-off distance, and the crosslinking limit is reached.
Next, the chemical crosslinking process of EPDM is briefly explained, based on the research of Wang et al. [5,19], to model the process with a third 5-ethylidene-2-norbornene (ENB) monomer.van Duin et al. [20] and Papanikolaou et al. [2] showed that peroxide free radicals mainly capture the tertiary hydrogen (C2x) of the main chain and the hydrogens at the C3 and C9 positions of ENB, similar to Papanikolaou et al.The positions are shown schematically in Figure 2a.In EPDM, the main chain-to-side-chain ratio is approximately 9:1.Zachary et al. [23] showed that the reaction rate of C3 in the monomer was 90%, while that of C9 was only 10%.The C2x atom of the main chain and the C3 atom of the side chain will be where the crosslinking occurs.The chemical crosslinking with the C9 atom was omitted because its reaction rate is quite low compared to the previous ones.
pairs are exhausted within the pre-assigned cut-off distance, and the crosslinking limit reached.
Next, the chemical crosslinking process of EPDM is briefly explained, based on th research of Wang et al. [5,19], to model the process with a third 5-ethylidene-2-norbornen (ENB) monomer.van Duin et al. [20] and Papanikolaou et al. [2] showed that peroxid free radicals mainly capture the tertiary hydrogen (C2x) of the main chain and the hydr gens at the C3 and C9 positions of ENB, similar to Papanikolaou et al.The positions a shown schematically in Figure 2a.In EPDM, the main chain-to-side-chain ratio is appro imately 9:1.Zachary et al. [23] showed that the reaction rate of C3 in the monomer wa 90%, while that of C9 was only 10%.The C2x atom of the main chain and the C3 atom the side chain will be where the crosslinking occurs.The chemical crosslinking with th C9 atom was omitted because its reaction rate is quite low compared to the previous one Physical crosslinking occurs through weak interactions.These interactions are due electrostatic and van der Waals forces that influence the gelation and elasticity of the po ymer.Physical crosslinking in elastomers largely determines the mechanical properties the elastomer and provides stability to the material.Figure 2b  After the initial models are sufficiently equilibrated to simulate the EPDM curin process, as described in the previous paragraph, curing is added periodically by impl menting a stepwise algorithm or multi-step subroutine as described below: Step 1: all di tances between reactive atoms were calculated and crosslinks were created between pai that were within a predetermined cutoff distance.Step 2: continuously, the covalent term of the formed crosslinks were created and the hydrogen atoms were removed from th newly formed pairs.Step 3: the new bonds were turned on in a stepwise manner by alte nating the bond equilibrium distances and force coefficients to avoid atoms missing fro the simulation box due to large atomic forces.Step 4: next, an NPT simulation was pe formed and the simulation box was further equilibrated before the distances between th reactive atoms were checked again.Step 5: the procedure was repeated until there we no other reactive pairs within the pre-assigned cutoff distance.
Subsequently the EPDM simulation box was deformed.Uniaxial stress simulation in the unit cell are performed by increasing the length along the loading direction (strai Physical crosslinking occurs through weak interactions.These interactions are due to electrostatic and van der Waals forces that influence the gelation and elasticity of the polymer.Physical crosslinking in elastomers largely determines the mechanical properties of the elastomer and provides stability to the material.Figure 2b schematically shows chemical crosslinking and physical crosslinking.
After the initial models are sufficiently equilibrated to simulate the EPDM curing process, as described in the previous paragraph, curing is added periodically by implementing a stepwise algorithm or multi-step subroutine as described below: Step 1: all distances between reactive atoms were calculated and crosslinks were created between pairs that were within a predetermined cutoff distance.Step 2: continuously, the covalent terms of the formed crosslinks were created and the hydrogen atoms were removed from the newly formed pairs.Step 3: the new bonds were turned on in a stepwise manner by alternating the bond equilibrium distances and force coefficients to avoid atoms missing from the simulation box due to large atomic forces.Step 4: next, an NPT simulation was performed and the simulation box was further equilibrated before the distances between the reactive atoms were checked again.Step 5: the procedure was repeated until there were no other reactive pairs within the pre-assigned cutoff distance.
Subsequently the EPDM simulation box was deformed.Uniaxial stress simulations in the unit cell are performed by increasing the length along the loading direction (strain) at each MD step.A value of 1 × 10 8 s −1 was assigned for the strain rate, which is a typical value in MD simulations due to the small-time scale set and is impossible to implement in a real physical system according to experimental standards.The volume of the simulation box during the deformation process is constant, the EPDM is considered incompressible, with a Poisson's ratio equal to 0.5, and the pressure in the transverse directions is considered constant, using a barostat to regulate it.The strain limit is determined when it reaches a value of 200%.The stress-strain curves were obtained for the different temperature sets (ranging from 428 K to 478 K with 10 K steps).
Viral expression is a widely used method in atomistic calculations for the intramolecular stress tensor [24][25][26], given in Equation (8), and this method provides a measure of the mechanical interactions between molecules in a continuous form: where index i runs over all particles 1 through N; m i , v i , and f ij denote the mass, velocity and force acting on particle i; and V denotes the volume of the system.To calculate the i − th contribution of forces f i acting on the system from position r i , the notation r ij is used, which indicates the position and the particular force that contributes to the total sum of forces f ij on i.

Phenomenological Models
The literature cites several mathematical models to describe elasticity in rubbery materials from a phenomenological approach.These models have constants that are estimated from experimental data by nonlinear regression.The stress-strain curves can be fitted using the semi-empirical formulation of the Mooney-Rivlin model.
The Mooney-Rivlin model can have 2, 3, 5, and 9 parameters to describe the behavior of rubber-like materials [27,28].The potential strain energy function, W, of an isotropic material, is usually formulated in terms of three invariants of the strain relations.The principal invariants, I 1 , I 2 , and I 3 , of the right Cauchy-Green deformation tensor are considered.The deformation energy potential functions are assumed to have polynomial or reduced polynomial forms.The polynomial model for an incompressible elastomer can be expressed in general terms as follows in Equation ( 9) [29]: where W is the work of deformation per unit volume or the elastically stored free energy density, C ij are the empirically determined material constants and I 1 and I 2 are the first and second Green's deformation invariants, respectively.For incompressible rubber materials, the expansion of the Mooney-Rivlin model for two parameters can be expressed as follows [30].
Green's deformation invariants are related to strain ratios (SR) in the three main directions by: For incompressible materials, it is considered that where λ 1 , λ 2 and λ 3 are the strain ratios in the three principal directions, considering the incompressibility ratio for uniaxial stress [31]: The crosslinking and entanglement effects of elastomeric polymers are represented by the model parameters C 1 and C 1 , respectively.
In rubbery polymers, elasticity is determined by the crosslink density resulting from the vulcanization process of the elastomer.For isotropic incompressible elastomers, the elastic stored free energy density, induced by the strain-energy relationship, is given by the specific phenomenological model of the Mooney-Rivlin Equation (11).
Given these relationships, Equation (11) becomes: The Kirchhoff stress tensor and the Green's strain tensor are used to determine the general state of the tension with the following equation: where t ij is the Kirchhoff stress tensor and γ ij is Green strain tensor.
For an incompressible material subjected to uniaxial tension, the stress σ is determined by deriving Equation (12) in terms of the strain ratio: Taking dW/dλ = σ into account and rearranging it, we finally have: To determine the parameter C 1 , Equation ( 16) is used.This equation includes the crosslink density N c of the vulcanized EPDM [7,13].
In the physics of entangled networks, there is currently no sufficiently rigorous and constructive formalism that adequately describes entanglement for networks of macroscopic size.Several researchers have proposed that C 2 takes on a value of zero because there is no formal equation for determining the parameter C 2 [32].
Given the above considerations, the Mooney-Rivlin model can accurately describe the mechanical behaviour of EPDM in components subjected to less than 200% deformation.

Molecular Dynamics Simulation
The crosslink density (N c ) and crosslink rate .N c for each of the test temperatures, Figures 3a and 3b, respectively, were obtained from the MD simulation of the EPDM composite.These results are plotted on a normalized time scale (τ); they show that the density and rate are strongly temperature-dependent.The ratio of the number of crosslinks formed to the probability of the total possible bond formation defines the degree of vulcanization and will be a function of the active positions in the system, electrostatic interactions, van Ver Waals forces and entanglement chains.The change in the slope of the crosslinking rate curves marks the gel point at the transition of the amorphous material (see Figure 4a).This transition at the gel point is due to increased crosslinking, which limits molecular movement through the polymer network.At the molecular level, the vulcanization process is characterized by chain polymerization and crosslinking, which cause significant changes in the mechanical properties, as shown macroscopically.Figure 4b shows the evolution of the density and rate of crosslink formation as a function of temperature.At higher temperatures, the polymer network becomes saturated with crosslinks, which affects the elasticity of the elastomer.At the molecular level, the vulcanization process is characterized by the formation of chains and crosslinking, which causes substantial changes in the mechanical properties that are observed macroscopically.At higher temperatures, the polymeric network becomes saturated with crosslinks, which will have an impact on the resulting properties of the elastomer.An important parameter is the mean squared displacement (MSD).Figure 5a shows a clear trend of the MSD of EPDM decreasing with the increasing temperature.This can be attributed to the combined effect of two key factors: the restrictions on chain movement due to crosslinking and the strong interactions at the bond interface.The change in the slope of the crosslinking rate curves marks the gel point at the transition of the amorphous material (see Figure 4a).This transition at the gel point is due to increased crosslinking, which limits molecular movement through the polymer network.At the molecular level, the vulcanization process is characterized by chain polymerization and crosslinking, which cause significant changes in the mechanical properties, as shown macroscopically.Figure 4b shows the evolution of the density and rate of crosslink formation as a function of temperature.At higher temperatures, the polymer network becomes saturated with crosslinks, which affects the elasticity of the elastomer.The change in the slope of the crosslinking rate curves marks the gel point at the transition of the amorphous material (see Figure 4a).This transition at the gel point is due to increased crosslinking, which limits molecular movement through the polymer network.At the molecular level, the vulcanization process is characterized by chain polymerization and crosslinking, which cause significant changes in the mechanical properties, as shown macroscopically.Figure 4b shows the evolution of the density and rate of crosslink formation as a function of temperature.At higher temperatures, the polymer network becomes saturated with crosslinks, which affects the elasticity of the elastomer.At the molecular level, the vulcanization process is characterized by the formation of chains and crosslinking, which causes substantial changes in the mechanical properties that are observed macroscopically.At higher temperatures, the polymeric network becomes saturated with crosslinks, which will have an impact on the resulting properties of the elastomer.An important parameter is the mean squared displacement (MSD).Figure 5a shows a clear trend of the MSD of EPDM decreasing with the increasing temperature.This can be attributed to the combined effect of two key factors: the restrictions on chain movement due to crosslinking and the strong interactions at the bond interface.At the molecular level, the vulcanization process is characterized by the formation of chains and crosslinking, which causes substantial changes in the mechanical properties that are observed macroscopically.At higher temperatures, the polymeric network becomes saturated with crosslinks, which will have an impact on the resulting properties of the elastomer.An important parameter is the mean squared displacement (MSD).Figure 5a shows a clear trend of the MSD of EPDM decreasing with the increasing temperature.This can be attributed to the combined effect of two key factors: the restrictions on chain movement due to crosslinking and the strong interactions at the bond interface.As crosslinking increases and electrostatic interactions occur, the chains begin to in tertwine, resulting in the stiffening of the geometrically constrained network.The move ment (diffusion) of the polymer chains is restricted, and the confinement leads to the co existence of an amorphous crosslinked structure forming a gel state.Figure 6a illustrate the evolution of the radius of gyration in crosslinking, which is associated with the stiff ness of the resulting system that manifests itself with the progressive increase in elasticity At the temperature of 428 K, there is greater flexibility in the main chain, which i expected to favor the formation of crosslinks.At 478 K, an agglomeration of side chains i formed.There is no significant conformational change in the main chain due to the low availability of free volume.
The radial distribution function denoted by g(r) defines the probability of finding certain distance between two atoms of different polymeric chains where there is a physica connection by electrostatic forces.Figure 6b represents the radial distribution functio where the electrostatic equilibrium distance (physical bonds) is 1.1 Angstrom, represente in the figure with the dotted red line.The value of the radial function increases with tem perature.
In the stress-strain ratio curves obtained by stretching in the MD simulation (Figur 7a,b), an increase in stiffness is observed concerning the increase in the tested temperature The increase in stiffness is related to the ability of the polymer chains to slide with eac As crosslinking increases and electrostatic interactions occur, the chains begin to intertwine, resulting in the stiffening of the geometrically constrained network.The movement (diffusion) of the polymer chains is restricted, and the confinement leads to the coexistence of an amorphous crosslinked structure forming a gel state.Figure 6a illustrates the evolution of the radius of gyration in crosslinking, which is associated with the stiffness of the resulting system that manifests itself with the progressive increase in elasticity.As crosslinking increases and electrostatic interactions occur, the chains begin to tertwine, resulting in the stiffening of the geometrically constrained network.The mo ment (diffusion) of the polymer chains is restricted, and the confinement leads to the existence of an amorphous crosslinked structure forming a gel state.Figure 6a illustra the evolution of the radius of gyration in crosslinking, which is associated with the st ness of the resulting system that manifests itself with the progressive increase in elastici At the temperature of 428 K, there is greater flexibility in the main chain, which expected to favor the formation of crosslinks.At 478 K, an agglomeration of side chains formed.There is no significant conformational change in the main chain due to the lo availability of free volume.
The radial distribution function denoted by g(r) defines the probability of findin certain distance between two atoms of different polymeric chains where there is a physi connection by electrostatic forces.Figure 6b represents the radial distribution functi where the electrostatic equilibrium distance (physical bonds) is 1.1 Angstrom, represent in the figure with the dotted red line.The value of the radial function increases with te perature.
In the stress-strain ratio curves obtained by stretching in the MD simulation (Figu 7a,b), an increase in stiffness is observed concerning the increase in the tested temperatu The increase in stiffness is related to the ability of the polymer chains to slide with ea At the temperature of 428 K, there is greater flexibility in the main chain, which is expected to favor the formation of crosslinks.At 478 K, an agglomeration of side chains is formed.There is no significant conformational change in the main chain due to the low availability of free volume.
The radial distribution function denoted by g(r) defines the probability of finding a certain distance between two atoms of different polymeric chains where there is a physical connection by electrostatic forces.Figure 6b  In the stress-strain ratio curves obtained by stretching in the MD simulation (Figure 7a,b), an increase in stiffness is observed concerning the increase in the tested temperature.The increase in stiffness is related to the ability of the polymer chains to slide with each other (the higher the crosslinking density, the lower the mobility) and the rearrangement of the polymer chains in the available free volume fraction.The free volume is the space not occupied by the polymer chains within the unit c and decreases with the increasing crosslink density and chain entanglement.
The lower the free volume fraction (FFV), the more difficult it is for the chains rearrange, limiting their mobility, which can be seen macroscopically as an increase stiffness.
Understanding stress behavior in terms of free volume implies that for conform tional changes to occur in the backbone, there must be space available to which the m lecular segment can move, implying a deformation that can be observed macroscopica Due to the conformational rearrangement induced by uniaxial deformation, wh the minimum value of the FFV is reached, the chains can no longer slide, thus restrict the molecular motion, which is observed as an increase in the slope of the S-s plot.
The concept of free volume explains the dependence of the mobility or conform tional rearrangement of molecular chains on temperature and plays an important role the mechanics of polymers.
The FFV can be calculated using Equation (17) where  is the free volume and is the volume occupied by the molecular chains in the unit cell.The free volume is the space not occupied by the polymer chains within the unit cell and decreases with the increasing crosslink density and chain entanglement.
The lower the free volume fraction (FFV), the more difficult it is for the chains to rearrange, limiting their mobility, which can be seen macroscopically as an increase in stiffness.
Understanding stress behavior in terms of free volume implies that for conformational changes to occur in the backbone, there must be space available to which the molecular segment can move, implying a deformation that can be observed macroscopically.
Due to the conformational rearrangement induced by uniaxial deformation, when the minimum value of the FFV is reached, the chains can no longer slide, thus restricting the molecular motion, which is observed as an increase in the slope of the S-s plot.
The concept of free volume explains the dependence of the mobility or conformational rearrangement of molecular chains on temperature and plays an important role in the mechanics of polymers.
The FFV can be calculated using Equation ( 17) where V f is the free volume and V o is the volume occupied by the molecular chains in the unit cell.
Molecular transport occurs when a chain moves into voids larger than the critical volume.The cooperative motion of neighboring atoms creates the voids by redistributing the free volume.The free volume provides further insight into the mechanical properties of EPDM, and Figure 8a-f shows the influence of the temperature on the distribution of the occupied volume as it increases.

𝑉 + 𝑉
Molecular transport occurs when a chain moves into voids larger than the critical volume.The cooperative motion of neighboring atoms creates the voids by redistributing the free volume.The free volume provides further insight into the mechanical properties of EPDM, and Figure 8a-f shows the influence of the temperature on the distribution of the occupied volume as it increases.Table 1 shows the calculated values of the free volume fraction (Equation ( 17)) for the different simulation temperatures.

Phenomenological Model
Given the results obtained from the MD simulations, the Mooney-Rivlin (MR) model, Equation (15), was employed to determine the material constants.
Equation ( 18) minimizes the total squared error between the data obtained from the MD simulation and the data fitted by the MR model, which give  pairs of data  ,  ,  = 1, … , , to determine  and  .

𝑒 = 𝜎 (𝜆 ) − 𝜎 (18)
Figure 9a shows the stress obtained from the MD simulation and the stress calculated using the material parameters  and  fitted with the MR model.Figure 9b   Table 1 shows the calculated values of the free volume fraction (Equation ( 17)) for the different simulation temperatures.

Phenomenological Model
Given the results obtained from the MD simulations, the Mooney-Rivlin (MR) model, Equation (15), was employed to determine the material constants.
Equation ( 18) minimizes the total squared error between the data obtained from the MD simulation and the data fitted by the MR model, which give N pairs of data {λ i , σ i }, i = 1, . . ., N, to determine C 1 and C 2 .
Figure 9a shows the stress obtained from the MD simulation and the stress calculated using the material parameters C 1 and C 2 fitted with the MR model.Figure 9b shows the C 1 and C 2 trends from the Mooney-Rivlin model fitted to the MD simulation data.It is widely documented that the C 1 parameter depends on the crosslinking density.
Figure 9a shows the stress obtained from the MD simulation and the stress calculated using the material parameters  and  fitted with the MR model.Figure 9b  To calculate the temperature-dependent crosslink density, N c (T), the MD simulation results (Table 1) are fitted to an order-two polynomial model, Equation (19).The proposed model is limited to the range (350 K-480 K), where the curing rate is greater than zero and the temperature is below the point where physical curing (due to electrostatic interactions) begins to degrade.
where N c (T) is the temperature-dependent crosslink density, a 0 , a 1 and a 2 are the parame- ters of the tested material and T is the temperature within the range to evaluate the crosslink density.
In this research, a modified empirical model, Equation (19), was proposed to calculate the C 1 (T) Equation ( 20) parameter of the Mooney-Rivlin model, which takes into account the temperature-dependent variation of the crosslinking density N c (T).
where the units of C 1 are Pa, R in MPa•cm 3 /mol•K, crosslink density N c (T) in mol/cm 3 and the absolute temperature (T) in Kelvin.Figure 10a shows a maximum limit of the possible crosslink formation density and, therefore, the existence of the maximum tensile stress before the reversibility or breakage of the physical crosslinks occurs due to the increase in temperature.Figure 10b   To calculate the temperature-dependent crosslink density,  (), the MD simulation results (Table 1) are fitted to an order-two polynomial model, Equation (19).The proposed model is limited to the range (350 K-480 K), where the curing rate is greater than zero and the temperature is below the point where physical curing (due to electrostatic interactions) begins to degrade.

𝑁 (𝑇) = 𝑎 + 𝑎 𝑇 + 𝑎 𝑇
where  () is the temperature-dependent crosslink density,  ,  and  are the parameters of the tested material and  is the temperature within the range to evaluate the crosslink density.
In this research, a modified empirical model, Equation (19), was proposed to calculate the  () Equation ( 20) parameter of the Mooney-Rivlin model, which takes into account the temperature-dependent variation of the crosslinking density  ().
where the units of  are  , R in MPa • cm mol • K ⁄ , crosslink density  () in mol cm ⁄ and the absolute temperature (T) in Kelvin.Figure 10a shows a maximum limit of the possible crosslink formation density and, therefore, the existence of the maximum tensile stress before the reversibility or breakage of the physical crosslinks occurs due to the increase in temperature.Figure 10b  Several investigations report that the  parameter is related to entanglements, weak interactions and the volume fraction occupied by the polymer networks.
The mechanical strength of polymers is significantly affected by the conformational reorganization of the molecular chains as a function of temperature.
An empirical model proposed in this research to determine the  parameter in Equation ( 21) includes the volume fraction occupied by the increase in the formation of Several investigations report that the C 2 parameter is related to entanglements, weak interactions and the volume fraction occupied by the polymer networks.
The mechanical strength of polymers is significantly affected by the conformational reorganization of the molecular chains as a function of temperature.
An empirical model proposed in this research to determine the C 2 parameter in Equation ( 21) includes the volume fraction occupied by the increase in the formation of physical crosslinks.At the macroscopic level, a change in the gradient of the stress-strain ratio occurs.
V O (T) is the fraction of the occupied volume as a function of temperature described by an Arrhenius-type model, where A is the pre-exponential factor and E is the specific activation energy, and this term comes from the phenomenological model of curing (Isayev-Deng) obtained by fitting the conversion fraction curve [33][34][35], a term widely used in the study of rubber vulcanization.
where V O (T) is the occupied volume fraction as a function of temperature and S 0 is a fitting parameter in MPa. Figure 11a graphically shows that as the temperature increases, the formation of crosslinks increases proportionally because the expansion of the linked chains within the unit cell tends to occupy the available volume (FFV). () is the fraction of the occupied volume as a function of temperature described by an Arrhenius-type model, where  is the pre-exponential factor and  is the specific activation energy, and this term comes from the phenomenological model of curing (Isayev-Deng) obtained by fitting the conversion fraction curve [33][34][35], a term widely used in the study of rubber vulcanization.
where  () is the occupied volume fraction as a function of temperature and  is a fitting parameter in MPa. Figure 11a graphically shows that as the temperature increases, the formation of crosslinks increases proportionally because the expansion of the linked chains within the unit cell tends to occupy the available volume ().
Figure 11b shows the parameter  versus temperature.This plot shows the contrast of the parameter  obtained from the fit with the MR model and the proposed model as a function of the occupied volume fraction.2 (Equations ( 19), ( 21) and ( 22)).Lengger et al. [11] correlated molecular-level crosslinking and chain entanglement as the main factors causing significant changes in the macroscopically observed mechanical properties of the polymer, i.e., increased stiffness.They also found that for parameters  and  , there is a strong correlation with the reduction in the available volume, and these show exponential saturation when they reach approximately 88.6% of their maximum values.19), ( 21) and ( 22)).
Lengger et al. [11] correlated molecular-level crosslinking and chain entanglement as the main factors causing significant changes in the macroscopically observed mechanical properties of the polymer, i.e., increased stiffness.They also found that for parameters C 1 and C 2 , there is a strong correlation with the reduction in the available volume, and these show exponential saturation when they reach approximately 88.6% of their maximum values.In their research, Panyukov et al. [32] analyzed structures with highly overlapping and interconnected loops of a finite size, the conformation of combined chains in polymer networks and their influence on the resulting mechanical properties.In this investigation, in agreement with previous reports by Lengger and Panyukov and their respective collaborators, it is found that the parameters of the Mooney-Rivlin model depend strongly on the volume fraction of the polymer occupied by the crosslinking and entanglement of the chains.Table 3 shows the RM model parameters obtained from Equations ( 15), ( 20) and (22).

Crosslinking Density, Mechanical Stress and Their Relation
The temperature increase is proportional to the amount of energy required to achieve the deformation, according to the results of molecular dynamics simulations.This effect is manifested macroscopically as an increase in stiffness or the Young's modulus (see Figure 7b).According to Saleesung et al. [36], Xie et al. [37], Papanikolaou et al. [2] and Bandyopadhyay et al. [38], with increasing crosslinking and decreasing free volume, the Young's modulus and thus the tensile strength of the elastomer increases.Wang et al. [19] reported that the effect of crosslinking between the backbones or physical crosslinking increases the stiffness modulus and decreases the free volume fraction.Gomez et al. [4] found that at higher temperatures, the interactions of the molecules are weakened by the accelerated movement of the radius of gyration in the segments, causing the physical crosslinks to break at some point.
The mechanical properties of the polymer are related to moving molecular chains and varying free volumes.When mechanical stress is applied to EPDMs, the free volume distribution changes as the polymer chains reorganize due to stretching.This process continues until a steady state is reached.The increase in the crosslinking density is responsible for the expansion of the molecular chains within the unit cell, which implies a decrease in the free volume fraction, resulting in the increased stiffness and tensile strength of the EPDM elastomer.
For deformations greater than 300%, the MD simulation data no longer fit the Mooney-Rivlin model.By increasing the deformation up to 800%, significant changes in the S-s ratio can be observed in Figure 12a, such as the tensile strength, which increases with the increasing crosslink density but tends to stabilize at a given crosslink density.The mechanical strength of the cured rubber is also affected by the entanglement of the molecular chains.Rubbers with higher physical crosslinking (lower FFV) tend to have the highest tensile strength values, while materials with higher FFV have the lowest tensile strength values (see Figure 12a,b).Having said the above, Wang et al. [5] and Dijkhuis et al. [39] agree that free energy and the amount of physical crosslinking are related to changes in the stiffness of a material.In plot 11b, for the S-s curves, in the strain ratio range above 600%, it is observed at  = 428 K , the tensile strength does not stabilize.At  = 468 K , the tensile stren reaches a stable maximum value, while at  = 478 K, the behavior can be interpreted fracture or failure of the material.
Physical tests have been carried out to verify the effect of the temperature on hardness obtained during the injection molding process in the door grommet (see The IRHD scale is defined as follows: if the Young's modulus is zero, the IRHD ha ness will be 0 and a material with an infinite Young's modulus will correspond to an IR hardness of 100.Rubber hardness ranges from 30 to 90 IRHD, with 45 to 65 conside adequate for door grommet applications and 65 to 85 for bumper applications.Tab shows the hardness test results of the manufactured samples.In plot 11b, for the S-s curves, in the strain ratio range above 600%, it is observed that at T = 428 K, the tensile strength does not stabilize.At T = 468 K, the tensile strength reaches a stable maximum value, while at T = 478 K, the behavior can be interpreted as a fracture or failure of the material. Physical tests have been carried out to verify the effect of the temperature on the hardness obtained during the injection molding process in the door grommet (see Figure 13a).Hardness tests were performed on the cured part by ASTM 1415-88 [40] using a PCE-DX-AS Rubber and Elastomer Hardness Tester (Shore A hardness 0-100+) and an International Rubber Hardness Degree Test (IRHD), with a measurement scale between 0 and 100, whose value is non-linear and related to the Young's modulus (modulus of longitudinal elasticity) of the material.The International Hardness test is based on measurements of the penetration of a rigid ball into the rubber specimen under specified conditions.Figure 13b shows the areas of measurement.In plot 11b, for the S-s curves, in the strain ratio range above 600%, it is observed that at  = 428 K , the tensile strength does not stabilize.At  = 468 K , the tensile strength reaches a stable maximum value, while at  = 478 K, the behavior can be interpreted as a fracture or failure of the material.
Physical tests have been carried out to verify the effect of the temperature on the hardness obtained during the injection molding process in the door grommet (see Figure 13a).Hardness tests were performed on the cured part by ASTM 1415-88 [40]  The IRHD scale is defined as follows: if the Young's modulus is zero, the IRHD hardness will be 0 and a material with an infinite Young's modulus will correspond to an IRHD hardness of 100.Rubber hardness ranges from 30 to 90 IRHD, with 45 to 65 considered adequate for door grommet applications and 65 to 85 for bumper applications.Table 4 shows the hardness test results of the manufactured samples.The IRHD scale is defined as follows: if the Young's modulus is zero, the IRHD hardness will be 0 and a material with an infinite Young's modulus will correspond to an IRHD hardness of 100.Rubber hardness ranges from 30 to 90 IRHD, with 45 to 65 considered adequate for door grommet applications and 65 to 85 for bumper applications.Table 4 shows the hardness test results of the manufactured samples.
For soft materials, the hardness is between 10 and 35 IRHD and for tires, between 85 and 100 IRHD.The results of the hardness tests (ASTM 1415-88 standard) on the Door Grommet component show that the hardness, and, therefore, the Young's modulus, increases with the increasing temperature, which is consistent with the results obtained from the stress-strain ratio curves of the molecular dynamics simulation (see Figure 7a,b), where the change in the slope with the temperature can be seen, as well as with the data adjusted to the Mooney-Rivlin model (see Figure 9a).The results of the MR model can be correctly extrapolated to the Door Grommet applications since the deformation range does not exceed 200% in use.
The results of this research indicate that for elastomeric polymers, there is a temperature (T u ) at which a threshold crosslink density is reached.In an injection molding process below the threshold temperature, an under-cure phenomenon will occur, while in a process with a temperature above T , an over-cure phenomenon will occur.
If under-cure occurs, the consistency will be soft and may not return to its original dimensions when deformed over 600%.When over-cure occurs, it will crack when subjected to deformations above 600%.When curing occurs at temperatures near T u and the strain ratio exceeds 600%, the tensile strength reaches a stable maximum.
The mechanical behavior is reasonably consistent with the Mooney-Rivlin model for strain ratios of less than 300% in the temperature range tested.

Conclusions
The results presented in this research provide new molecular insights into the crosslinking mechanisms of elastomeric polymers and their influence on mechanical behavior, which could facilitate the design of rubbery products in applications such as door grommets.In combination with recent publications on tensile strength behavior using the Mooney-Rivlin model for rubbery polymers, they demonstrate the ability of MD simulations to elucidate complex phenomena associated with deformation and crosslinking behaviors in polymers.
The time-scale used in the molecular dynamics simulations might seem to be a limitation for the superposition of results; the phenomenological model presented in this research to predict the stress-strain behavior depends on the crosslink density, which can be understood macroscopically as vulcanization, so the time-scale is not an important issue.In this case, the molecular dynamics results allow us to perform multi-scale simulations.
This research provides evidence that the temperature during the curing process plays a preponderant role in the formation of crosslinks, resulting in a complex mixture of crosslinked structures (chemical and physical), as well as the movement of the backbone of the EPMD polymer trying to form clusters of crosslinks between the main chains, generating a significant conformational change, and a decrease in the free volume which translates into greater stiffness at higher temperatures.
The proposed model for determining Mooney-Rivlin parameters in elastomeric materials, which takes into account the effect of varying crosslink densities and volume fractions occupied by polymer chains as a function of temperature, provides a practical approach to analyze the impact on the mechanical performance of rubbery materials in applications such as door grommets and other types of seals.However, it is not limited to this, as any

Figure 1 .
Figure 1.(a) Random distribution of EPDM polymer chains in a simulation box (each color represent different backbone).(b) Procedure for molecular dynamics simulation.

Figure 1 .
Figure 1.(a) Random distribution of EPDM polymer chains in a simulation box (each color represent different backbone).(b) Procedure for molecular dynamics simulation.

Figure 2 .
Figure 2. (a) Scheme of physical crosslinking rubber with peroxide of ENB-EPDM (+ENB-EPD and H transfer or + macro-radical (EPDM)).(b) Schematic representation of the molecular structu of a cured elastomer.(c) Chemical structure of ethylene-propylene-diene monomer.

Figure 2 .
Figure 2. (a) Scheme of physical crosslinking rubber with peroxide of ENB-EPDM (+ENB-EPDM and H transfer or + macro-radical (EPDM)).(b) Schematic representation of the molecular structure of a cured elastomer.(c) Chemical structure of ethylene-propylene-diene monomer.

Figure 4 .
Figure 4. (a) Rate of crosslink density versus Crosslink density.(b) Variation of Crosslink Density and Crosslink Density Rate with Temperature.

Figure 4 .
Figure 4. (a) Rate of crosslink density versus Crosslink density.(b) Variation of Crosslink Density and Crosslink Density Rate with Temperature.

Figure 4 .
Figure 4. (a) Rate of crosslink density versus Crosslink density.(b) Variation of Crosslink Density and Crosslink Density Rate with Temperature.

Figure 6 .
Figure 6.(a) Radius of gyration evolution for the trajectories of the different temperatures.(b) Radial distribution function, (or pair correlation ().

Figure 6 .
Figure 6.(a) Radius of gyration evolution for the trajectories of the different temperatures.(b) R dial distribution function, (or pair correlation function) ().

Figure 6 .
Figure 6.(a) Radius of gyration evolution for the trajectories of the different temperatures.(b) Radial distribution function, (or pair correlation function) represents the radial distribution function where the electrostatic equilibrium distance (physical bonds) is 1.1 Angstrom, represented in the figure with the dotted red line.The value of the radial function increases with temperature.

Figure 7 .
Figure 7. (a) Stress-strain ratio behavior evolution curve.(b) Stress-strain ratio curve between m imum and minimum simulation temperatures.
when a chain moves into voids larger than the criti volume.The cooperative motion of neighboring atoms creates the voids by redistribut the free volume.The free volume provides further insight into the mechanical propert of EPDM, and Figure8a-f shows the influence of the temperature on the distribution the occupied volume as it increases.

Figure 7 .
Figure 7. (a) Stress-strain ratio behavior evolution curve.(b) Stress-strain ratio curve between maximum and minimum simulation temperatures.
Figure9ashows the stress obtained from the MD simulation and the stress calculated using the material parameters  and  fitted with the MR model.Figure9bshows the  and  trends from the Mooney-Rivlin model fitted to the MD simulation data.It is widely documented that the  parameter depends on the crosslinking density.

Figure 9 .
Figure9ashows the stress obtained from the MD simulation and the stress calculated using the material parameters  and  fitted with the MR model.Figure9bshows the  and  trends from the Mooney-Rivlin model fitted to the MD simulation data.It is widely documented that the  parameter depends on the crosslinking density.
graphically compares C 1 of the MR model and C 1 (T) of the proposed model.

Materials 2024 , 20 Figure 9 .
Figure 9. (a) Tensile strength curves obtained by Molecular Dynamics and curves fitted to the Mooney-Rivlin model.(b) Evolution of the parameters of the Mooney-Rivlin model as a function of temperature.

Figure 10 .
Figure 10.(a) Fitting of the proposed crosslink density model.(b) Contrasting the  parameter of the MR model with that calculated using the proposed exponential model (the blue dotted lines represent the limits of the studied range).

Figure 10 .
Figure 10.(a) Fitting of the proposed crosslink density model.(b) Contrasting the C 1 parameter of the MR model with that calculated using the proposed exponential model (the blue dotted lines represent the limits of the studied range).

Figure 11 .
Figure 11.(a) Volume fraction distribution versus temperature.(b) Contrasting the  parameter of the MR model with that calculated using the proposed exponential model.Niu et al. [3] proposed an exponential-type model to determine the parameters of rubber material  and  .They obtained good results with a mean deviation between 4.1% and 5.2% from the calculated results and combined with nonlinear autoregression (NAR).The degradation of crosslinking as a function of time is studied indirectly in the Niu et al. model; in the case of this investigation, the RM parameters are a function of temperature.The phenomenology is consistent because the empirical relationship between the two models focuses on studying elastomer crosslink evolution.With the model proposed in this research, a mean deviation of less than 3.5% was obtained.The values of the fitting parameters are shown in Table2(Equations (19), (21) and (22)).

Figure 11 .
Figure 11.(a) Volume fraction distribution versus temperature.(b) Contrasting the C 2 parameter of the MR model with that calculated using the proposed exponential model.

Figure
Figure 11b shows the parameter C 2 versus temperature.This plot shows the contrast of the parameter C 2 obtained from the fit with the MR model and the proposed model as a function of the occupied volume fraction.Niu et al. [3] proposed an exponential-type model to determine the parameters of rubber material C 1 and C 2 .They obtained good results with a mean deviation between 4.1% and 5.2% from the calculated results and combined with nonlinear autoregression (NAR).The degradation of crosslinking as a function of time is studied indirectly in the Niu et al. model; in the case of this investigation, the RM parameters are a function of temperature.The phenomenology is consistent because the empirical relationship between the two models focuses on studying elastomer crosslink evolution.With the model proposed in this research, a mean deviation of less than 3.5% was obtained.The values of the fitting parameters are shown inTable 2 (Equations (19), (21) and (22)).Lengger et al.[11] correlated molecular-level crosslinking and chain entanglement as the main factors causing significant changes in the macroscopically observed mechanical properties of the polymer, i.e., increased stiffness.They also found that for parameters C 1 and C 2 , there is a strong correlation with the reduction in the available volume, and these show exponential saturation when they reach approximately 88.6% of their maximum values.

Figure 13 .
Figure 13.(a) Injection molding process in the door grommet; (b) Hardness testing measurement areas.

Figure 13 .
Figure 13.(a) Injection molding process in the door grommet; (b) Hardness testing measurement areas.

Table 1 .
Crosslinking density, crosslinking rate and free volume fraction of EPDM compound at different temperatures.

Table 1 .
Crosslinking density, crosslinking rate and free volume fraction of EPDM compound at different temperatures.